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ABSTRACT 

We perform numerical simulations on the merger of multiple black holes (BHs) in 
primordial gas at early cosmic epochs. We consider two cases of BH mass: Mbh = 
30 Mq and Mbh = 10 4 M 0 . Attention is concentrated on the effect of the dynamical 
friction by gas in a host object. The simulations incorporate such general relativistic 
effects as the pericentre shift and gravitational wave emission. As a result, we find 
that multiple BHs are able to merge into one BH within 100 Myr in a wide range of 
BH density. The merger mechanism is revealed to be categorized into three types: gas- 
drag-driven merger (type A), interplay-driven merger (type B), and three-body-driven 
merger (type C). We find the relation between the merger mechanism and the ratio of 
the gas mass within the initial BH orbit (M gas ) to the total BH mass (Y Mbh)- Type 
A merger occurs if M gas > 10 5 Y Mbh, type B if M gas < 10 5 Y Mbe;, and type C if 
M gas <C 10 5 Y -Mbh- Supposing the gas and BH density based on the recent numerical 
simulations on first stars, all the BH remnants from first stars are likely to merge into 
one BH through the type B or C mechanism. Also, we find that multiple massive BHs 
(Mbh = 10 4 Mq) distributed over several parsec can merge into one BH through 
the type B mechanism, if the gas density is higher than 5 x 10 6 cm -3 . The present 
results imply that the BH merger may contribute significantly to the formation of 
supermassive BHs at high redshift epochs. 

Key words: dark ages, reionization, first stars - stars: black holes - galaxies: high- 
redshift - quasars: supermassive black holes - gravitational waves - methods: numer¬ 
ical 


1 INTRODUCTION 

In the last two decades, it has been demonstrated that 
massive galaxies harbor supermassive black holes (SMBHs) 
in the centres of their bulge components (Kormendy & 
Ho 2013, and references therein). Also, at redshifts higher 
than six, quasars are found that possess SMBHs with the 
mass higher than 10 9 Mq (Fan et al. 2001; Kurk et al., 
2007). The formation history of these SMBHs is among the 
most significant unsolved issues in astrophysics (Volonteri 
2012; Haiman 2013). Two recently discovered high-redshift 
quasars, ULASJ112010+641 with the mass of Mbh = 2 x 
10 9 M 0 at redshift z = 7.085 (Mortlock et al. 2011) and 
SDSS J01001+2802 with M BH = 1.2 x 10 10 M© at z = 6.30 
(Wu et al. 2015) have raised a serious problem for the for¬ 
mation of SMBHs. Possible building blocks for such high- 
redshift SMBHs are the remnants of first stars. The ini¬ 
tial mass function of first stars is thought to be more or 
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less top-heavy (Abel et al. 2000; Nakamura & Umemura 
2001; Bromm et al. 2002; Yoshida et al. 2006; Greif et al. 
2011; Susa et al. 2014; Hirano et al 2014). First stars of sev¬ 
eral tens M© undergo supernovae, leaving black holes (BHs) 
of few tens M© (Heger & Wooseley 2002). For SMBHs to 
grow from such first star remnants through mass accretion 
at z > 6, a super-Eddington accretion rate is requisite. If 
SMBHs grow continuously by mass accretion from BH rem¬ 
nants of ~ 20 Mq, the Eddington ratio (A) is required to 
be A = 1.4 for ULASJ112010+641, or A = 1.3 for SDSS 
J01001+2802. However, the continuous accretion is unlikely 
to be sustained due to feedbacks, and thus the average mass 
accretion rates should be lower than the Eddington rate 
(Milosavljevic 2009; Alvarez et al 2009). On the other hand, 
seed BHs may stem from supermasssive stars of 10 4 ~ 6 M© 
as a result of the direct collapse of primordial density fluc¬ 
tuations (Umemura, Loeb & Turner 1993; Bromm & Loeb 
2003; Inayoshi & Omukai 2012). These BHs are thought to 
be incorporated into a primordial galaxy of ~ 10 8 — 10 9 M© 
(Greene 2012). If an SMBH grows via gas accretion from 
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such a massive BH, the constraint on the accretion rate can 
be alleviated. 

Another possible pathway of SMBH formation is the 
merger of BHs. If the merger of multiple black holes pre¬ 
cedes the growth via gas accretion, the merged BH can be 
a seed of a supermassive black hole. So far, the merger of 
SMBHs in a galaxy has been argued extensively. As for a 
binary of SMBHs, Begelman et al. (1980) pointed out that 
the orbit of a binary SMBH cannot shrink below one par¬ 
sec due to the loss cone depletion (the depletion of stars 
on orbits that intersect the binary SMBH), which is often 
called the final parsec problem (e.g. Merritt et al. 2004). 
It is argued that if the host galaxy provides an aspherical 
potential, a binary SMBH may overcome the final parsec 
problem (Khan, Just & Merritt 2011, Khan et al. 2012; 
Khan et al. 2013). However, this solution of the problem 
is still under debate (Vasiliev, Antonini & Merritt 2014). If 
there are more than two SMBHs in a galaxy, the dynamical 
relaxation of SMBHs is significantly controlled by the grav¬ 
ity of SMBHs themselves, especially by three-body interac¬ 
tion. When a third MBH intrudes into a binary SMBH, one 
SMBH carries away angular momentum from the remain¬ 
ing two SMBHs, reducing the binary separation and even¬ 
tually inducing the merger of the binary (Iwasawa, Funato 
& Makino 2006). In the case of many SMBHs, the stellar dy¬ 
namical friction allows a binary MBH to interact frequently 
with other SMBHs, and then the decay of the binary orbits 
leads to the merger (Tanikawa & Umemura 2011, 2014). 

In a first-generation object formed at an early cosmic 
epoch, the dynamical friction by stars is unlikely to work 
effectively, since the initial mass function is top-heavy and 
most stars undergo supernovae. However, the dynamical fric¬ 
tion by gas could work. Recent radiation hydrodynamic sim¬ 
ulations on the formation of first stars show that several or 
more stars form in a primordial gas cloud with the den¬ 
sity of around 10 7 cm -3 and the extension of 1000 AU, 
where the gas fraction is 99% (Greif et al. 2011; Umemura 
et al. 2012; Susa 2013; Susa et al. 2014) . In this circum¬ 
stance, BH remnants of first stars are most likely subject to 
the dynamical friction by abundant gas. The gas dynamical 
friction has been considered as a mechanism that prompts 
the BH merger (Ostriker 1999, Tanaka & Haiman 2007). 
Hitherto, the merger processes by the gas dynamical fric¬ 
tion have been investigated in the case of two massive BHs 
(e.g. Escala 2004, Escala 2005). In this paper, we explore 
the merger of multiple BHs, supposing a first-generation ob¬ 
ject of ~ 10 s — 10 6 Mq or a gas-rich primordial galaxy of 
- 10 s - 10 9 M e . 

The paper is organized as follows: In section 2, we de¬ 
scribe the numerical method. In section 3, we show the nu¬ 
merical results. In section 4, we discuss the merger criterion 
through the gas friction. In section 5, we summarize the pa¬ 
per. 


2 METHOD OF NUMERICAL SIMULATIONS 

Here, we present the framework of numerical simulations. 


2.1 Equation of motion 

The equations of motion for BHs are given by 


d 2 n 
dt 2 



+ a PN,i? 


+ a DF,i T a po t,i, 


( 1 ) 


where r,; and r ; are respectively the positions of i -th BH and 
j -th BH, Abh is the number of BHs, G is the gravitational 
constant, rrij is the mass of j -th BH, apN.ij is the general 
relativistic acceleration of j-th BH on 1-th BH in the Post- 
Newtonian prescription, a|j^ i is the dynamical friction (DF) 
on 1-th BH by gas, and a po t,i is the acceleration on 1-th BH 
by gravitational potential of gas. 


2.2 Key parameters 

We initially set up multiple BHs of equal mass. In this paper, 
we consider two cases of the BH mass: one is 30 Mg BHs 
as first star remnants, which are born in a first-generation 
object of ~ 10 s — 10 6 Mq, and the other is 10 4 Mq BHs 
resulting from supermassive stars, which are incorporated 
into a primordial galaxy of ~ 10 8 —10 9 Mg. A key parameter 
in the present simulations is the BH density, Pbh, at the 
initial epoch. Recent radiation hydrodynamic simulations on 
first star formation have shown that several or more stars are 
born in a disk of ~ 1000 AU (Greif et al. 2011; Umemura et 
al. 2012; Susa 2013; Susa et al 2014). We change the typical 
extension of the BH distribution at the initial epoch, rt yp , 
to settle the BH density. In the case of 30 Mq mass BH, we 
alter rt yp from 0.01 pc to 1 pc. In the case of 10 4 Mg mass 
BH, we set rt yp from 0.1 pc to 10 pc. Then, the BH density 
is given by 


Pbh = 


477 r ty P 


Y, mi ' 


( 2 ) 


where rm is the mass of i-th BH. Another key parameter is 
the gas number density n gas . Simulations of first stars show 
that the density in a first object is ~ 10'~ 8 cm~ 3 in the 
star-formation epoch. Here, we consider a wider range of 
the gas density from 10 4 cm' 3 to 10 12 cm -3 to elucidate 
the dependence of the BH merger on the gas density. 


2.3 Gas drag and potential 

We give the gas dynamical friction and the gas potential by 
analytic solutions. We use the formula of the gas dynamical 
friction force given by Tanaka & Haiman (2009) for the mo¬ 
tion with Mi < M eq and Ostriker (1999) for Mi > M eq , 
where Mi is the Mach number of i-th BH and M eq is the 
Mach number where these two formulas give equal accel¬ 
eration. Here, we adopt M eq = 1-5 as Tanaka & Haiman 
(2009). Then, the acceleration of the gas dynamical friction 
( a DF,i) is s iven b y 

a DF,i = -MG 2 mimHn gaB (r)^xf(Mi) 

Vi 


( 3 ) 
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where mu is the mass of the hydrogen atom, n gas is the 
number density of gas, Vi is the velocity of i-th BH, and t 
is the elapsed time. The r m j n is the minimum scale of the 
dynamical friction on a BH, and we give r m i n as Gmi/v 2 . 
Here, Vit means the effective scale of gas medium, and we 
set an upper limit of Vit to 0.1 pc. When Vit < r m j n , we 
assume f(Mi) = 0. 

In this paper, we postulate uniform background gas to 
purely extract the dependence on the gas density, and give 
its density as a parameter. Then, the gravitational accelera¬ 
tion by gas (a potj i) and its time derivative (a po t,i) are given 
as 

4 

a po t,i — “TrGmHUgasrj, (5) 

4 

a po t,i — — 7rGmHngasVi. (6) 


2.4 Gas temperature 


The heating rate due to the gas dynamical friction is esti¬ 
mated as follows (Kim et al. 2005) : 

Abp = 7 -3xl0- 25 ergcm- 3 s- 1 ( I5; ^) 

M bh V > 

30 M© ) V 2 

T \~ 1/2 / WBH 
1000 I<) V 10 pc- 3 

where hbh is the number density of BHs, and the angular 
brackets denote the average over the Maxwellian distribu¬ 
tion, 





f(v) 


47tA^bH 2 —u 2 /(2 <t 2 ) 

-4i ' v r> 

(27rcr 2 ) 3 / 2 


( 8 ) 


where ay is the one-dimensional velocity dispersion. 

In a first object, the cooling is dominated by hydrogen 
molecules (H 2 ) at T k, 10 3 K (e.g. Omukai 2000). I 11 a pri¬ 
mordial galaxy, the temperature goes down by the cooling 
of neutral hydrogen (HI) around T ~ 10 4 K, and further re¬ 
duces by the H 2 cooling down to T « 10 3 K, if the gaseous 
metallicity is lower than a percent of solar abundance (Susa 
& Umemura 2004). The H 2 cooling rate through H — H 2 col¬ 
lision (Hollenbach & McKee 1989) and the HI cooling rate 
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(Thoul & Weinberg 1995) are respectively given by 


A H2 (10 3 K) = 2.8 x lo¬ 

ng 

x 


erg cm s 
/h 2 


)( 


10 4 cm -3 7 V 3 x 10 


(9) 


A H i(10 4 K) 


1.0 x 10~ 14 erg cm -3 s _1 


10 4 cm -3 


) 2 /i i 


( 10 ) 


where /h 2 is the fraction of H 2 molecules and /hi is the 
fraction of neutral hydrogen atoms. Here, /h 2 ^ 3x 10~ 4 
in the range of n gas > 10 4 cm" 3 (Palla, Salpeter & Stahler 
1983). 

We find that, in both cases of BH mass, the heating rate 
is lower than the cooling rate over the range of parameters 
in which the simulations are performed. Therefore, the gas 
temperature is expected to settle at T ~ 10 3 K. In this 
paper, we assume the gas temperature to be 1000 K. Then, 
the sound speed is given as C s = 3.709 [km/s]. 


2.5 Relativistic effects 


We incorporate the general relativistic effects according to 
the Post Newtonian prescription up to 2.5PN term (Kupi 
2006). The relativistic effects on i-th BH by j- th BH are 
given by 1PN (aiPN.ij) and 2PN (a 2 PN,ij) term correspond¬ 
ing to the pericentre shift, and 2.5PN term (a2.5pN.ij) cor¬ 
responding to the gravitational wave (GW) radiation. They 
are expressed by 


aiPN.ij 


a2PN,ij 


Grrij 

[»[- 

r?- 

+s ( 

Gnu 


+ (vi 

~ Vj 

Grrij 

[n[—: 

r 2 

r ij 
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V r i: 
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a2.5PN,ij 


where 
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rij = (r* - r j). 


(15) 


2.6 Merger condition 

We assume that two MBHs merge, when their separation is 
less than 100 times the sum of their Schwarzschild radii: 

|ri —rj | < 100 (r sch ,i + r sc h,j ), (16) 

where r sc h,i is the Schwarzschild radius of j-th BH given 
by 2 Gmi/c 2 with the speed of light c. At the final stage of 
merger, the binding energy of a binary BH is transformed 
to the energy of GW, retaining the mass of each BH. Hence, 
the BH mass after the merger is just the sum of two MBHs. 




i[10 6 yr] 


2.7 Numerical scheme 


We integrate the equations of motion using the fourth- 
order Hermite scheme with the shared time step (Makino 
& Aarseth 1992). To use the Hermite scheme, we calculate 
the time derivative of the acceleration by the Newtonian 
gravity, the gas gravitational potential, and the relativistic 
force. On the other hand, we treat the dynamical friction of 
gas to quadratic order, since the formula of the dynamical 
friction is not so accurate as that of gravity and also the 
fourth-order scheme is especially requisite during the grav¬ 
itational three-body interaction between a close BH binary 
and an intruding BH. Therefore, the time derivative of the 
gas dynamical friction is not included. Then, the shared time 
step in the Hermite scheme is given by 


At = mini 


UHer,! 11 a Her,z 

\ V \a. W lla (3) 

\ l cl Her,i I l c *’Her,i 


+ l a Her,i I” 


+ l a Her,i P 


(17) 


where aHer,i is the acceleration on z-th BH that is treated 
to the fourth-order, and given by 


^BH 

<3-Her,i = ^ ^ a GlTLj 


+ apN ,ij f + a po t,i- 


(18) 


(k) 

The a^ er i is the fc-th derivative of aHer.i, and y is the accu¬ 
racy parameter. We assume y = 0.003 in the present simu¬ 
lations. 

To avoid cancellation of significant digits when such tiny 
scales as 100 r sc h are resolved, we calculate the BHs evolu¬ 
tion in the coordinate where the origin is always set to the 
center of mass for the closest pair of BHs. This prescription 
is based on the simulations on the merger of multiple BHs 
(Tanikawa & Umemura 2011), and allows us to pursue ac¬ 
curately the orbit of the BHs until the merger condition is 
satisfied. 


Figure 1. The separation of the closest pair within all BHs as 
a function of time. The initial BH mass is Mbh = 30 Mg, the 
initial typical extension of the BH distribution is rt yp = 1.0 pc, 
and the gas density is n gas = 5 X 10 7 (top) or 5 X 10® (bottom) 
cm -3 . 


2.8 Setup of simulations 

We set up ten BHs as a fiducial case, and besides investigate 
the cases of two or three BHs to scrutinize the key physics 
of merger. The initial positions of BHs are given randomly 
in the x — y plane within r typ . Also, we give the velocity 
to each BH as the sum of a circular component and a ran¬ 
dom component. The circular velocity is given to balance 
against the gravity in the x — y plane. The random velocity 
is given according to a Gaussian distribution with the same 
dispersion as the circular velocity. The random velocity is 
given in the xyz space. In the case of two BHs, we give only 
circular velocity without initial eccentricity. This condition 
is desired to discriminate the BH merger purely by the gas 
friction. 

We calculate each run until 100 Myr, since the back¬ 
ground environments of the host objects are likely to change 
in 100 Myr. Also, we terminate the simulation, if all BHs 
merge into one BH. 


3 NUMERICAL RESULTS 

3.1 Dependence on gas density and BH density 

The results of ten BHs systems are summarized in Table [T| 
and Table [2] Table |T] is those for the case of Mbh = 30 Mg 
and Tabic [2] is for Mbh = 10 4 Mg. The column is the ini¬ 
tial extension of the BH distribution (rt yp ) and the corre¬ 
sponding BH density (pbh), while the row is the assumed 
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Table 1. The results with Mbh = 30 Mg and Abh = 10. 


r typ [pc] 

1.0 

0.464 

0.215 

0.1 

0.0464 

0.0215 

0.01 

Pbh[Mqpc- 3 ] 

7.2 x 10 1 

7.2 x 10 2 

7.2 x 10 3 

7.2 x 10 4 

7.2 x 10 s 

7.2 x 10 6 

7.2 x 10 7 


A m 

type 

A m 

type 

Am 

type 

Am 

type 

Am 

type 

Am 

type 

Am 

type 

n gas [cm -3 ] 

tfln[yr] 

tfin[yr] 

tfin[yr] 

tfin [yr] 

tfin[yr] 

tfin[yr] 

tfin[yr] 

10 12 

0 

5 

A 

9 

A 

9 

A 

9 

A 

9 

A 

9 

B 

1.0 x 10“ 

1.0 x 10° 

3.6 x 10' 

4.4 X 10° 

1.7 x 10° 

7.0 x 10° 

4.3 x 10° 

10 11 

2 

A 

8 A 

9 

A 

9 

A 

9 

A 

9 

A 

9 

B 

1.0 x 10° 

1.0 x 10° 

1.3 x 10' 

1.8 x 10° 

1.8 x 10° 

2.3 x 10“ 

2.8 x 10“ 

10 10 

5 

A 

9 

A 

9 

A 

9 

A 

9 

A 

9 

B 

9 

B 

1.0 x 10° 

4.5 x 10' 

5.1 x 10° 

5.9 x 10° 

7.4 x 10“ 

6.3 x 10“ 

1.7 x 10° 


5 

A 

9 

A 

9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

5 x 10 9 

1.0 x 10° 

2.9 x 10' 

3.1 x 10° 

4.5 x 10° 

1.6 x 10° 

3.4 x 10° 

7.7 x 10“ 


8 A 

9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

9 

B 

10 9 

1.0 x 10° 

2.4 X 10' 

3.5 x 10° 

2.6 x 10° 

4.3 x 10° 

3.5 x 10° 

3.0 x 10° 

5 x 10 s 

8 A 

9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

9 

B 

1.0 x 10“ 

1.2 X 10' 

1.3 x 10° 

6.5 x 10° 

4.3 x 10° 

5.5 x 10° 

5.1 x 10° 


9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

10 8 

3.0 x 10' 

5.1 x 10° 

4.0 x 10° 

5.5 x 10° 

3.6 x 10° 

4.2 x 10° 

1.2 x 10' 


9 

A 

9 

B 

9 B 

9 

B 

9 

B 

9 

B 

9 

B 

5 x 10 7 

4.5 x 10' 

3.7 x 10° 

2.2 x 10' 

3.2 x 10' 

1.3 x 10' 

4.7 x 10° 

3.6 x 10° 


9 

B 

9 

B 

9 B 

9 

B 

9 

B 

9 

B 

9 

B 

10 7 

3.8 x 10' 

2.3 x 10' 

1.7 x 10' 

3.3 x 10' 

1.8 x 10' 

2.9 x 10' 

1.7 x 10' 


9 

B 

9 

B 

9 B 

9 

B 

9 

B 

9 

B 

9 

B 

5 x 10® 

4.2 X 10' 

3.9 x 10' 

4.2 x 10' 

4.7 x 10' 

6.3 x 10' 

3.5 x 10' 

3.1 x 10' 


6 

B 

6 

B 

8 B 

6 

C 

8 C 

6 

C 

6 

C 

10® 

1.0 x 10° 

1.0 x 10° 

1.0 x 10 s 

1.0 x 10° 

1.0 x 10“ 

1.0 x 10° 

1.0 x 10° 


2 

C 

6 

C 

6 

C 

4 

C 

5 

C 

3 

C 

4 

C 

5 x 10 5 

1.0 x 10° 

1.0 x 10° 

1.0 x 10° 

1.0 x 10° 

1.0 x 10° 

1.0 x 10° 

1.0 x 10° 


0 

0 

1 

C 

0 

1 

C 

2 

C 

0 

10 5 

1.0 x 10° 

1.0 x 10“ 

1.0 x 10° 

1.0 x 10 s 

1.0 x 10° 

1.0 x 10° 

1.0 x 10“ 


0 

0 

0 

1 

C 

0 

0 

0 

10 4 

1.0 x 10“ 

1.0 x 10° 

1.0 x 10 s 

1.0 x 10° 

1.0 x 10“ 

1.0 x 10° 

1.0 x 10° 


gas density (n gas ) in the system. N m is the number of the 
merged BHs (nine means all BHs merged into one). The 
’type’ shows the merger mechanism, which is classified in 
the next section. Also, the termination time of simulations, 
tfin, is shown. 

In the case of Mbh = 30 Mq (Tabic [TJ , we find that 
if tigas is in the range from 5 x 10 6 to 10 s cm -3 , then all 
the BHs merge into one massive BH over six orders of BH 
density (pbh = 72 — 7.2 x 10' Mgpc -3 ). Also, if the gas 
density is higher than n gas = 5 x 10 5 cm -3 , the BH merger 
occurs once at least in 100 Myr. In the recent numerical 
simulations on first stars (Greif et al. 2011; Susa 2013; Susa 
et al. 2014), the density of BH remnants is expected to be 
Pbh ~ 10' Mq pc -3 and the gas density is to be in the range 
of n gas « 10' — 10 s cm -3 . Hence, the present numerical 
results imply that all the BH remnants from first stars are 
likely to merge into one BH. 

In the case of Mbh = 10 4 Mq (Table [2]), we find that if 


n gas is in the range from 5x 10® to 10 9 cm -3 , then all the BHs 
merge into one massive BH over five orders of BH density 
(p B H = 2.4 x 10 2 — 2.4 x 10 7 Mq pc -3 ). Interestingly, in the 
gas denser than 5 x 10® cm -3 , multiple BHs can merge into 
one, even if the typical separation is greater than several pc. 
This result implies that the BH merger is able to contribute 
significantly to the formation of SMBHs at high redshifts. 
Similar to the case of Mbh = 30 Mq, if the gas density is 
higher than n gas = 5 x 10 5 cm -3 , the BH merger occurs once 
at least. 

In both cases of BH mass, if the gas density is very high 
and simultaneously the BH density is very low (upper left in 
Table [T] and Table H, no merger occurs. This is due to the 
fact that the deep gravitational potential by gas enhances 
the circular velocities of BHs and therefore a self-gravitating 
system of BHs is hardly able to form. 

As shown in terms of the types of merger in Table [T] and 
Table [2] the merger mechanism changes systematically with 
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Table 2. The results with Mbh = 10 4 A/q and JVbh = 10. 


rtyp [PC] 

10.0 

4.64 

2.15 

1.0 

0.464 

0.215 

0.1 

Pbh[MqPc- 3 ] 

2.4 x 10 1 

2.4 x 10 2 

2.4 x 10 3 

2.4 x 10 4 

2.4 x 10 s 

2.4 x 10 6 

2.4 x 10 7 


N m 

type 

N m 

type 

Am 

type 

A m 

type 

N m 

type 

N m 

type 

N m 

type 

rigas [cm -3 ] 

tfin[yr] 

tfln[yr] 

tfin[yr] 

*fin [yr] 

tfin[yr] 

tfln[yr] 

tfin[yr] 

10 12 

0 

4 

A 

7 

A 

9 

A 

9 

A 

9 

A 

9 

A 

1.0 x 10 s 

xl0 lu 

xl0 lu 

1.3 x 10' 

1.5 x 10° 

1.7 x 10° 

2.2 X 10" 

10 11 

0 

5 

A 

9 

A 

9 

A 

9 

A 

9 

A 

9 

B 

1.0 x 10 s 

1.0 x 10° 

4.6 x 10' 

5.0 x 10° 

5.8 x 10° 

7.2 x 10" 

1.6 x 10" 

10 10 

4 

A 

8 A 

9 

A 

9 

A 

9 

A 

9 

A 

9 

B 

1.0 x 10° 

1.0 x 10 s 

1.7 x 10' 

2.0 x 10° 

2.5 x 10° 

9.3 x 10" 

4.1 x 10" 


5 

A 

9 

A 

9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

10 9 

1.0 x 10° 

5.8 x 10' 

7.0 x 10° 

9.6 x 10 a 

7.6 x 10° 

9.0 x 10° 

2.4 X 10° 


5 

A 

9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

9 

B 

5 x 10 8 

1.0 x 10° 

5.0 x 10' 

1.0 x 10' 

1.5 x 10° 

9.5 x 10° 

7.6 x 10° 

6.7 x 10° 


7 

A 

9 

A 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

10 8 

1.0 x 10° 

2.5 x 10' 

3.7 x 10° 

8.3 x 10° 

4.9 x 10° 

5.6 x 10° 

7.9 x 10° 


9 

A 

9 

A 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

5 x 10 7 

7.3 x 10' 

2.0 x 10' 

1.3 x 10' 

7.1 x 10° 

1.2 x 10' 

3.5 x 10° 

9.4 x 10° 


9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

10 7 

7.7 x 10' 

3.5 x 10' 

9.3 x 10' 

6.0 x 10' 

3.3 x 10' 

3.6 x 10' 

4.1 x 10' 


8 B 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

9 

B 

5 x 10 6 

1.0 x 10 s 

8.3 x 10' 

3.9 x 10' 

8.5 x 10' 

4.6 x 10' 

4.8 x 10' 

7.9 x 10' 


5 

B 

7 

B 

4 

B 

9 

B 

5 

C 

5 

C 

6 

C 

10 6 

1.0 x 10° 

1.0 x 10° 

1.0 x 10 s 

6.5 x 10' 

1.0 x 10° 

1.0 x 10° 

1.0 x 10“ 


3 

B 

3 

C 

3 

C 

3 

C 

6 

C 

4 

C 

4 

C 

5 x 10 5 

1.0 x 10° 

1.0 x 10° 

1.0 x 10 s 

1.0 x 10“ 

1.0 x 10° 

1.0 x 10° 

1.0 x 10“ 


0 

0 

0 

0 

0 

1 

C 

0 

10 5 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10“ 

1.0 x 10 s 

1.0 x 10° 

1.0 x 10“ 


0 

0 

1 

C 

0 

1 

C 

1 

C 

1 

C 

5 x 10 4 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10° 

1.0 x 10° 

1.0 x 10“ 


0 

0 

0 

0 

0 

0 

1 

C 

10 4 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10 s 

1.0 x 10“ 

1.0 x 10 s 

1.0 x 10“ 

1.0 x 10“ 


the gas density (n gaa ) and the BH density (pbh). Details of 
the merger mechanism are described in the following. 

3.2 Types of merger mechanism 

Here, we scrutinize the merger mechanism and classify the 
mechanism. When a BH binary merges due to the dynam¬ 
ical friction by gas, the separation of the same pair of BHs 
monotonically shrinks, while the three-body interaction be¬ 
tween a close BH binary and an intruding BH often replaces 
one BH in the binary by the intruder and therefore the sep¬ 
aration of the closest pair changes violently. In Figs.[T][5l we 
show the separation of the closest pair within all BHs as a 
function of time, where the colors of lines change at every 
event of the BH merger. 

Fig.[T]shows the case of low-mass BHs (Mbh = 30 Mg). 
The initial extension of the BH distribution is rt yp = 1.0 pc. 
The upper panel corresponds to higher gas density, and 


the lower panel to lower gas density. As seen in the upper 
panel, the apocenter and pericenter distances of the clos¬ 
est pair decay smoothly, eventually resulting in the merger. 
Such smooth decay indicates that the gas dynamical friction 
drives the merger. On the hand, in the lower panel of Fig. 
□ the orbit is highly eccentric, since the apocenter and peri¬ 
center distances are largely different. The distances some¬ 
times oscillate violently. Such strong variance is thought to 
be caused by the three-body interaction, since the simul¬ 
taneous interaction with a 4th or 5th BH seems quite rare 
from a viewpoint of probability. Actually, we have confirmed 
that the contribution of a 4th or 5th BH is little. Thus, we 
call such interaction the three-body interaction. After the 
apocenter decays to less than ~ 10 -8 pc, the pair separa¬ 
tion lessens rapidly due to the effect of the GW emission, 
and eventually the BH binary merges. 

We classify the merger mechanism according to the 
manner of decay just before the GW promotes the merger. 
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Figure 2. Same as Figure [T] but for rt yp = O.lpc, and n gaa 
5 X 10® (top), 10 s (middle), or 5 X 10 s (bottom) cm -3 . 


Here, the merger mechanism is categorized into three types: 
gas-drag-driven merger (type A), three-body-driven merger 
(type C), and interplay-driven merger (type B). In type A, 
the dynamical friction by gas effectively decays the orbit and 
then the GW drives the merger. The type A is seen in the 
cases with higher gas density and lower BH density, as ex¬ 
pected. The examples of type A are shown in the top panels 
of Figs. 1 1 121 and HI and Fig. [3] In both types B and C, the 
strong disturbance of the orbit is induced by the three-body 
interaction during the first merger. But in type B, after the 
first few mergers by three-body interaction, the orbit decays 
slowly for long time through the gas dynamical friction. The 
examples of type B are seen in the middle panels of Figs. [2] 
and [4] and Fig. [5] In these results, after the mergers driven 
by three-body interaction, the orbit evolution starts from a 
larger separation than initial typical separation (rt yp ). Such 
increase of separation is due to the slingshot mechanism 
when an intruder interacts with the binary. This causes a 



Figure 3. Same as Figure [T] but for Mbh = 10 4 MQ,rt yp = 
10.0 pc, and ra gaa = 5 X 10 7 cm -3 . 


negative effect on the merger timescale. In other words, the 
three-body interaction may play a positive and negative role 
for the merger. 0 In type C, the strong disturbance of the 
orbit continues until the final merger. The three-body in¬ 
teraction of BHs solely transfers the angular momentum, 
eventually allowing the merger via the GW radiation. Type 
C is seen in the cases of high BH density. The examples of 
type C are shown in the bottom panels of Figs. 0and[|] 

As shown in Tables 1 and 2, in the cases of lower BH 
density (higher n yP ), type A is seen in a broad range of 
gas density. In the cases of higher BH density, type B is 
dominant in higher gas density environments. On the other 
hand, type C is common in lower gas density. But, if the 
gas density is very low, no merger occurs even in the same 
density of BHs. This implies that the gas dynamical friction 
plays a non-negligible role even for type C. 

3.3 BH number dependence 

To clarify the physics of BH merger, we perform reference 
simulations with the different number of BHs. In Fig. [6] 
we compare the resultant merger time in ten-BH systems 
(tAr=io) to that in three-BH (tjv =3 ) or two-BH systems 
{tN= 2 ) • Each panel corresponds to a different set of param¬ 
eters. Red, pink, and blue open symbols respectively rep¬ 
resent the types A, B, and C mergers in ten-BH systems. 
Brown- and green-filled symbols show the averaged merger 
time in the runs of three- and two- BH systems, respectively. 
Note that the averaged merger time (tfl n ) has the variance 
of about a factor of two that comes from the initial setup of 
random number. 

In the cases of two BHs, we give only circular velocity 
initially. Since there is no three-body interaction, the merger 
of two-BH systems is always driven by the gas dynamical 
friction. On the other hand, in the systems of three BHs, 
three-body interaction as well as the gas dynamical friction 

1 In the present study, we assume the uniform gas. But, in a real 
galaxy, the gas distribution is likely to become diffuse in outer 
regions. If the velocity of the kicked BH is high enough, the BH 
cannot fall back but may escape. 
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Figure 6. Merger time in ten-BH systems (t.y—io) as a function of gas number density n gas . The left panels are the cases with low-mass 
BHs (Mbh = 30 Mg), and right panels are those with high-mass BHs (Mbh = 10 4 Mq). The top, middle and bottom panels show the 
results for lower, intermediate, and higher BH density, respectively. Red, pink, and blue open symbols represent the results of the types 
A, B, and C mergers, respectively. Green- and brown-filled symbols are the results of two-BH (tiv= 2 ) and three-BH (by— 3 ) systems, 
respectively. 


can work to induce the merger. We can see that the merger 
timescales in type C accord with those in three-BH systems. 
In other words, no merger occurs in two-BH systems in the 
parameters of type C. Thus, the classification of type C in 
which the merger is driven by the three-body interaction is 
justified. 


On the other hand, in type A and type B, both two- 
BH and three-BH systems can merge. In the parameters 
of type A, the merger timescales in ten-BH systems are 
systematically shorter than those in three-BH systems. In 
the ten-BH systems, more BHs suffer from the gas dynami¬ 
cal friction and the angular momenta of BHs are extracted 
simultaneously. Therefore, the ten-BH systems can merge 
in shorter time on average. As shown in top left and top 
right panels, two-BH or three-BH systems cannot merge in 
higher gas density environments, whereas ten BHs systems 
can merge. These results imply that the dynamical friction 
on ten BHs enhances the three-body interaction, eventually 
allowing the merger. In the bottom right panel, the aver¬ 
aged merger time in ten-BH systems is slightly longer than 
that in three-BH systems, in the cases with lower gas den¬ 
sity (n gas < 10 8 cm -3 ). This is expected by the negative 
feedback as discussed in section 3.2. 


3.4 Effect of recoil 

In the above, we have not considered the effects of grav¬ 
itational wave recoil. The recoil velocity ranges from sev¬ 
eral ten km s _1 to several thousand km s _1 (Kesden et al. 
2010). Tanikawa and Umemura (2014) have shown that if 
the recoil velocity is lower than the virial velocity in the 
system, multiple BHs can eventually merge. In the present 
simulations, the virial velocity is higher than a few hun¬ 
dred km s -1 for n gas bs 10 8 cm -3 and M gas ^ 10 6 Mq in 
the case of Mbh = 30 Mq, and for n gas 2s 10® cm -3 and 
M gas ^ 10 7 Mq in the case of Mbh = 10 4 Mq. Actually, 
we have imposed random velocities on BHs at a level of the 
virial velocity, and revealed the conditions under which all 
BHs can merge. Hence, the merger condition is expected not 
to change, if the recoil velocity is lower than the virial ve¬ 
locity. To confirm this, we have simulated the several cases 
in which the recoil velocity is incorporated. As a result, we 
have found that if the recoil velocity is lower than the virial 
(escape) velocity, then the averaged merger time and the 
number of merged BHs are not changed. If the recoil ve¬ 
locity is higher than the escape velocity, then the merged 
BH escapes from the system. The actual magnitude of re¬ 
coil velocity is determined by the alignment of spins of two 
MBHs (Schnittman & Buonanno 2007; Kesden et al. 2010). 
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Figure 7. Same as Fig. [6] but the analytic estimates of gas dynamical friction timescale, / [>{,- , are plotted by green lines, instead of the 
results of two-BH and three-BH systems. 


The spins tend to be aligned through the gas accretion onto 
BHs. Assuming the level of the recoil velocity for BHs with 
aligned spins, we have performed further simulations and 
found that if the escape velocity of the system is higher 
than ~ 140km/s, then the BH merger proceeds in a similar 
way to the simulations without the recoil. 


4 CRITERION FOR GAS DRAG 

In Fig. [7] we compare the averaged merger time in ten BHs 
systems (tN=io) to the analytic estimate of gas dynamical 
friction timescale (Idf). We assess the friction timescale as¬ 
suming a binary BH, similar to those by Begelman et al. 
(1980) and Matsubayashi et al. (2007). Based on a test sim¬ 
ulation for a BH binary, we postulate that the eccentricity 
of the binary becomes very high, when the BH gravity dom¬ 


inates the gas potential. Then, Idf is estimated as 


toF = 

r l °’ {r) dr 


J 100 r sch r 

tDF(y) = 

^circ {fT ) 

maxlag^r')} 


„,3 

tDF(rtyp) — 

L/ circ 


... r ^ r < r typ 


(47rTn. H 7i gas ) 1/2 r^ Yp 

3 3 / 2 G 1 / 2 Mbh 

47rG 1 / 2 m H ngas7y ! (p 


(M gaSj b 3> Mbh) 
(M gaS! b « M B h) 


(19) 

( 20 ) 
( 21 ) 

( 22 ) 


Here, M gaH! b is the gas mass within the binary orbit, and 
Vcirc is the circular velocity given as i> c i rc = [G(2M gaSj b + 
MbH)/ 2r] 1 ^ 2 . In this estimate, we use the maximum value 
of a^p in the range of r ^ r' ^ r typ- Green curves in 
Fig- [3 present the estimated timescale of gas dynamical fric¬ 
tion. According to equation J22J, Idf is proportional to n ga s 
in the limit of high gas density, whereas £df is inversely 
proportional to n gas in the limit of low gas density. Thus, 
a turning point (the minimum) emerges in each curve of 
toF due to this change of dependence. As seen in Fig. [7] 
the turning point is in accordance with the transition from 
type A (gas-drag-driven merger) to type B (interplay-driven 
merger), and also the trend of merger time (t N=10 ) is in 
agreement with the analytic estimate of gas dynamical fric¬ 
tion timescale (Idf). The timescale becomes longer in the 
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Figure 4. Same as Figure □3 but for Mbh = 10 4 MQjT-typ = 
1.0 pc, and n gas = 10 9 (top), 10 8 (middle), or 5 x 10* (bottom) 
cm -3 . 


higher gas density due to the increase in the deepness of 
gravitational potential of gas. The longer timescale in the 
lower gas density stems from the strong gravitational poten¬ 
tial of other BHs. In addition, we can recognize the difference 
that the analytic estimate is systematically longer than the 
merger time in the simulations. One reason for the system¬ 
atic difference comes from dismissing the GW as well as the 
three-body interaction effect in the analytic estimate. An¬ 
other reason is concerns the fact that the analytic timescale 
is based on the approximate values of initial positions in 
type B and type C. 

In Fig. 0 we see the dependence of the merger timescale 
on the BH density (or r typ ). With the increasing BH density 
(decreasing r t yP ), the right side of the power law shifts lower. 
This shift stems from the decrease in the gravitational po¬ 
tential of gas within rt yp . This tendency is also seen in the 
numerical results of ten-BH systems. Thus, the rough ten¬ 
dency of the merger time in the type A region is explained 
by the analytic timescale. Therefore, we regard foF as an 



Figure 5. Same as Figure [T] but for Mbh = 10 4 Mg,rt yp = 
0.1 pc, and n gas = 10 8 cm^ 3 . 


appropriate estimate for the merger solely by the gas dy¬ 
namical friction. 

As described above, type A is always seen in higher gas 
density than the turning point of the double power law of 
Idf- These results imply that the turning point provides the 
critical gas density, above which the merger is driven by the 
gas dynamical friction. The turning point is determined by 
the ratio of the gas mass within rt yp , M gas , to the total BH 
mass, £M BH . In the analytic estimate, the turning point 
appears around M gas ~ 1.4 x 10 5 Mbh as the average of 
the each panel in Fig. [7] As for the numerical results, we 
can assess the turning points in the top and middle panels, 
since the turning points are not clear in the bottom panels. 
In the cases of 30 Mq BHs, the gas density at the turning 
points is 10 7 cm -3 in upper left panel or 10 10 cm~ 3 in middle 
left panels. Then, the ratio of M gas to ^ Mbh at the turning 
points is 1.4 x 10 5 in both cases. In the cases of 10 4 Mq BHs, 
it is 10' cm -3 in upper right panel, or 10 9 cm' 3 in middle 
right panels. Hence, the M gas -to-^ Mbh ratio is 4.2 x 10 5 or 
4.2 x 10 4 , respectively. These results show that the critical 
gas density above which type A is dominant is given by the 
condition of M gas / ^ Mbh ~ 10 5 . Based on the simulations 
on the formation of first objects (Greif et al. 2011; Susa et al. 
2014), M gas is likely to be lower than 10 s X) Mbh - Thus, type 
B and type C are expected to be more common mechanism 
of BH merger. 


5 CONCLUSIONS 

In this paper, we have investigated the merger of multiple 
BHs systems in an early cosmic epoch, incorporating the dy¬ 
namical friction by gas. For the purpose, we have performed 
highly accurate numerical simulations taking into account 
such general relativistic effects as the pericentre shift and 
gravitational wave emission. Consequently, we have found 
the following: 

1. Multiple BHs are able to merge into one BH within 
100 Myr in a wide range of parameters. In the case of 
Mbh = 30 Mq, if n gas is in the range from 5 x 10 6 
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to 10 8 cm~ 3 , then all the BHs can merge for the BH 
density of pbh = 72 — 7.2 x 10 7 Mgpc~ 3 . In the case of 
Mbh = 10 4 Mq, if n g as is in the range from 5 x 10 6 to 
10 9 end 3 , then all the BHs can merge for the BH density 
of pbh = 2.4 x 10 2 — 2.4 x 10' Mg pc -3 . Commonly, if the 
gas density is higher than n gas = 5 x 10 5 cm -3 , the BH 
merger occurs once at least in 100 Myr. 

2. The merger mechanism is classified into three types: 
gas-drag-driven merger (type A), three-body-driven merger 
(type C), and interplay-driven merger (type B). A criterion 
between type B and type C is almost concordant with the 
results of two- or three-BH systems. 

3. Based on the argument of merger timescales, type A 
merger has turned out to occur, if the gas mass within the 
initial BH orbit (M gas ) is higher than 10 5 J^.Mbh, where 
^2 Mbh is the total mass of BHs. 

4. Supposing the gas and BH density based on the recent 
numerical simulations on first stars (Greif et al. 2011; 
Umemura et al. 2012, Susa 2013; Susa et al. 2014), all the 
BH remnants from first stars are likely to merge into one 
BH through the type B or C mechanism. 

5. As for a primordial galaxy possessing massive BHs with 
Mbh = 10 4 Mq, it is suggested that in environments of gas 
density higher than 5 x 10 6 cm -3 , multiple BHs can merge 
into one through the type B mechanism, even if the typical 
separation is greater than several pc. This result implies 
that the BH merger is able to contribute significantly to the 
formation of SMBHs at high redshifts. 

A point we have not considered in this paper is the 
effect of mass accretion onto BHs. During the motion of 
BHs in gas, some gas can fall onto BHs through the Bondi- 
Hoyle-Littleton accretion. This effect may significantly affect 
the merger processes of BHs. This issue will be studied in 
a forthcoming paper. Besides, in a primordial galaxy, the 
stellar dynamical friction may work cooperatively with the 
gas friction to extract the angular momenta. Such synergism 
is also an issue to be explored in the future. 
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